FAST AND OBLIVIOUS CONVOLUTION QUADRATURE 
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Abstract. We give an algorithm to compute N steps of a convolution quadrature approximation 
to a continuous temporal convolution using only 0{N logN) multiplications and O(logAf) active 
memory. The method does not require evaluations of the convolution kernel, but instead 0{logN) 
evaluations of its Laplace transform, which is assumed sectorial. The algorithm can be used for the 
stable numerical solution with quasi-optimal complexity of linear and nonlinear integral and integro- 
differential equations of convolution type. In a numerical example we apply it to solve a subdiffusion 
equation with transparent boundary conditions. 
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tion, anomalous diffusion 
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1. Introduction. In this paper we give a fast and memory-saving algorithm for 
computing the approximation of a continuous convolution (possibly matrix x vector) 

t 

f{t~r)g{T)dr , < t < T, (1.1) 



by a convolution quadrature with a step size h > 0, 

n 

Y,^n-j gUh), n^l,...,N, (1.2) 

j=Q 

where the convolution quadrature weights ujn are determined from their generating 
power series as (see ^3 El El) 

f:^.c=F{'M). (1.3) 

n=0 

Here F{s) is the Laplace transform of the (possibly matrix-valued) convolution ker- 
nel fit), and (5(C) = 1 - C or d{C) = {1 - + - C? for the methods based 
on the first or second-order backward difference formula, respectively. We will also 
consider a similar approximation based on implicit Runge-Kutta formulas such as 
the Radau II A methods ^j. Attractive features of such convolution quadratures are 
that they work well for singular kernels f{t), for kernels with multiple time scales, and 
in situations where only the Laplace transform F{s) but not the convolution kernel 
fit) is known analytically. Perhaps most importantly, they enjoy excellent stability 
properties when used for the discretization of integral equations or integro-differential 
equations of convolution type, in a way often strikingly opposed to discretizations 
with more straightforward quadrature formulas (see references in \V2\). 
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The direct way to compute ()1.2|) is to first compute and store the (possibly matrix- 
valued) weights ujq, . . . ,wjv, which can be done accurately with 0{N) evaluations of 
the Laplace transform F{s) |TI], and then to compute the discrete convolution. Done 
naively, this requires 0{N'^) multiplications (possibly matrix x vector) and 0{N) 
active memory for the values g{jh) and for the weights. Using FFT, the number 
of multiplications can be reduced to 0{N logA^), and to 0(A^ (log A^)^) in the case 
of integral equations where the values of g{t) are not known beforehand, but where 
g{nh) is computed only in the nth time step However, that approach does not 
reduce the number of i^-evaluations and the memory requirements. 

Here we give an algorithm, also applicable in the case of linear and nonlinear 
integral equations, which computes (|1.2|l in a way that requires 

• 0{N \ogN) multiplications, 

• 0(log A^) evaluations of the Laplace transform F{s), and 

• 0(log A'') active memory. 

The history g{jh) for j = 0, . . . , is forgotten in this algorithm, and only logarithmi- 
cally few linear combinations of the g-values are kept in memory. These are obtained 
by solving numerically, with step size initial value problems of the form y' — Xy + g 
with complex A. The weights a;„ {n < N) are not computed explicitly, except the 
first few, e.g., the first 10 weights. 

The algorithm presented here uses ideas of the fast convolution algorithm of [T^. 
which instead of 1)1. 2|l makes a different approximation to the continuous convolution. 
The stability properties of the second-order method of ^3] in integro-differential equa- 
tions such as those of Section 5 are, however, extremely difficult to analyze (cf. also 
[ITp and remain entirely unclear for higher-order extensions. Here we show how 
the convolution quadratures (|1.2|l with all their known favorable properties can be 
implemented in an equally fast and memory-saving way. 

Following the error analysis of |Z1IH| we give exponentially convergent error bounds 
for the contour integral approximations that are employed in this algorithm. They 
ensure that the constants hidden in the O-symbols of the above work estimates depend 
only logarithmically on the error tolerance for these contour integral approximations. 

We assume a sectorial Laplace transform F{s): 

F{s) is analytic in a sector | arg(s — c)| < tt — with ip < ^tt, and there 
\F{s)\ < M Is]-" for some real M and > 0. 
The inverse Laplace transform is then given by 

fit) = 7^ [ e'^F{X)dX, t>0, (1.5) 

2711 Jp 

with F a contour in the sector of analyticity, going to infinity with an acute angle to 
the negative real half-axis and oriented with increasing imaginary part. The function 
f{t) is analytic in t > and satisfies 

|/(i)l <Ct''-ie^*, t>0, (1.6) 

and is therefore locally integrable. (The absolute values on the left-hand sides of 
the bounds (|1.4I) and (|1.6() are to be interpreted as matrix norms for matrix-valued 
convolution kernels.) 

In Section 2 we review convolution quadrature based on multistep and Runge- 
Kutta methods. We give a contour integral representation of the convolution quadra- 
ture weights whose discretization along hyperbolas or Talbot contours is discussed in 
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Section 3. The fast and oblivious convolution algorithm is formulated in Section 4. 
Finally, in Section 5 we give the results of numerical experiments with integral and 
integro-differential equations originating from regular and anomalous diffusion prob- 
lems. 

2. Convolution quadrature. In this section we review briefly convolution 
quadrature and give a contour integral representation of the convolution quadrature 
weights on which the fast algorithm of this paper is based. 

2.1. Convolution quadrature based on multistep methods. We consider 
the convolution quadrature l|1.2|l with weights (|1.3f) . By ()1.4|) and Cauchy's integral 
formula we have, with a contour F as in H1.5|) . 

Hence, with e„(z) defined by 

oo 

{5{C)-z)-' = Y.e.n{z)C, (2.1) 

we have the integral formula 

Wn = — en{h\)F{\)dX, (2.2) 

which can be viewed as the discrete analog of H1.5(l . For the backward Euler dis- 
cretization 5(C) = 1 — C we note the explicit formula 

e„(z) = (l-z)-"-\ (2.3) 

which is of the form e„(z) = g(z)r(z)" with r{z) = j-^ and q{z) — jz^- 

For the second-order BDF method, where d{() = X]fc=i 'h^^ ~ with p = 2, we 
obtain from a partial fraction decomposition of ((5(C) — z)^^ that 

e„(z) ^ -=i= f (2 - \/TT2^)-"-i - (2 + yTT2i)-"-i) , (2.4) 
Vl + 2z \ / 

which is of the form e„(z) — qi(z)ri(z)" + q2{z)r2{z)" . Connoisseurs of Cardano's 
formulas find analogous formulas to 12.4|l also for the BDF methods of orders 3 and 4. 

2.2. Convolution quadrature based on Runge-Kutta methods. We con- 
sider an implicit Runge-Kutta method with coefhcicnts aij, bj, for i, j = 1, . . . ,r7i. 
We denote the Runge-Kutta matrix by O = [o-ij), the row vector of the weights by 
6"^ = (6j), and the stability function by 

r{z) = 1 + zfe'^(/- za)-^ll, 

where U = (1, . . . , 1)^. We assume that all eigenvalues of the Runge-Kutta matrix d 
have positive real part and, for simplicity, that the method is A-stable and the row 
vector of the weights equals the last line of the Runge-Kutta matrix. 



hj = am] for j 1, . . . ,m, 
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and correspondingly c„i — 1. These conditions are in particular satisfied by the 
Radau IIA family of Runge-Kutta methods From such a Runge-Kutta method, a 
convolution quadrature is constructed as follows Let 

A(C) = ((3+^116^)"' (2.5) 

and define weight matrices Wn by 

f;H^„C"=^^(^). (2.6) 

Let w„ = (Wj^j, . . . , w™) denote the last row of Wn- Then an approximation to the 
convolution integral Hl.l|l at time i^+i ~ (n ^ l)/i is given by 

n ra n 

i=o i=l j=o 

with the column vector — (g{tj + Cih)y^^. For a Runge-Kutta method of classical 
order p and stage order q, this approximation is known to be convergent of the order 
min(p, q + 1 + f) with v of (|1.4|) . 

With the row vector e„(z) = (e^(z), . . . , e™(2:)) defined as the last row of the 
m X m matrix En{z) given by 

oo 

(A(C)-z/™)-i = 5]i?„(z)C", (2.8) 



n=0 



we obtain an integral formula like in (|2.2() . 

/" 

ujn^— en{h\) ^ F{X) dX. (2.9) 
27ri Jr 

For n > 0, e„{z) is given as 

e„(z) = r(z)Xz) (2.10) 
with the row vector q{z) = b^{I — zC2)~^; cf. Lemma 2.4 in [T3'. We note that 

n 

yi%^hJ2en^,{hX)g, (2.11) 

is the Runge-Kutta approximation at time tn+i of the linear initial value problem 

y' = Xy + g{t), y(0) = 0. (2.12) 
The convolution quadrature H2.7|l is thus interpreted as 

Un+i = 77- I vL+i dX ; 

see Proposition 2.1 in |13j . 



3. Approximation of the contour integrals. The fast convolution algorithm 
will be based on discretizing the integrals in H2.2|l and (|2.9|) along suitable complex 
contours. This approximation is discussed in the present section. 

3.1. Quadrature on Talbot contours and hyperbolas. The fast algorithm 
approximates the quadrature weights uJn by linear combinations of the exponential 
approximations en{hX), locally on a sequence of fast-growing time intervals nh € If. 



[B^-^h,2B%), 



(3.1) 



where the base _B > 1 is an integer. For example, B — \Q was found a good choice in 
our numerical experiments. The approximation on Ii results from applying the trape- 
zoidal rule to a parametrization of the contour integral for the convolution quadrature 
weights, 



27ri 



/ en{h\)®F{\)d\^h e„(/iAi^'^)(^F(Af ), nh e h, (3.2) 



k=-K 



with an appropriately chosen complex contour Yg. The number of quadrature points 
on r^, 2K -|- 1, is chosen independent of It is much smaller than what would be 
required for a uniform approximation of the contour integral on the whole interval 
[0, T]. Only a few of the first convolution quadrature weights, for n < No (e.g., 
A'o = 10), are approximated differently, using the trapezoidal rule discretization of 
the integral over a circle as discussed in j^J : 



= last row of / ^ 
./iChp 



n < Nn. 



(3.3) 



The numerical integration in H2.2|l or (|2.9() is done by applying the trapezoidal rule 
with equidistant steps to a parameterization of a hyperbola |H1 or a Talbot contour |18l 
The Talbot contour is given by 




Fig. 3.1. Talbot contour (left) and hyperbola (right). 



(-vr, TT)^r -.e^ 7(61) =a + ^i{e cot{e) + ivO) (3.4) 

where the parameters /i, v and a are such that the singularities of F{s) lie to the left 
of the contour and that the singularities of e„(ft,s) lie to the right of the contour. See 
left part of Figure IXTI for = 0. The parameter will depend on I via the right 
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end-point of J^, which yields a Talbot contour Ti depending on the approximation 
interval Ii. The weights and quadrature points in H3.2|l are given by (omitting i in 
the notation) 

"^fc = ~ r.^7-/i i\ I'i^k) , Afe=7(6'fc) with 6k 



2{K +1) ' ^ ' ' ' K+l 

Alternatively, the hyperbola is given by 

r : 61 7(61) = - sin(a + i6l)) (3.5) 

where the parameters /i and a are such that the singularities of F{s) lie to the left 
of the contour. See the right part of Figure ITTI for a = 7r/2 — 1/2. The weights and 
quadrature points in H3.2|) are given by (omitting I in the notation) 

"^fc = TT l'{(^k) , Afe = 7(6*^) with 9k ^ kr , 
Zn 

where r is a step length parameter. 

3.2. Numerical experiments. In view of the examples of Section|31we present 
here numerical experiments with 

fit) ^ for which F{s) = s-^/"^. 

V7rt 




10° 10^ 10^ 10^ 10^ 10° 10' 10^ 10° io' 



Fig. 3.2. Talbot quadrature errors versus time for K = 15, B = 5 and K = 10, B = 10 for 
different Integrators. (Implicit Euler, BDF(2), RadauIIA(3) and RadauIIA(5) in clockwise order 
starting from the upper left corner) 

The error is calculated with respect to a reference solution, obtained for a dis- 
cretization of the contour integral with a large number of integration points. For the 
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Fig. 3.3. Hyperbola quadrature errors versus time for K = 15, B = 5 and K = 10, B = 10 for 
different Integrators. (Implicit Euler, BDF(2), RadauIIA(3) and RadauIIA(5) in clockwise order 
starting from the upper left corner) 

Radau II A methods of order 3 and 5, where the LOn are row vectors of dimension 2 
and 3, respectively, we plot the error of the last entry. 

Using the Tabot contours, the following choices of parameters were found to give 
good results. A relative accuracy of about 10~^ on the interval h for t>2 with right 
end-point Ti is obtained with B — K = 10, ji — 8/Ti, v = 0.6. For a relative 
approximation error of 10~^, take i? = 5, if = 15, and the other parameters as before, 
cf. Fig. 13.21 For n > 20 there is no substantial difference between the different Runge- 
Kutta methods. Since the approximations to the first few convolution quadrature 
weights are poor, they will not be used in the algorithm. 

Using the hyperbola contours, a relative accuracy of about 10~^ on the interval Ii 
for ^ > 2 with right end-point is obtained with B — 10, K — 10, a — 1, jj. — 
and r = 0.64. For a relative approximation error of 3 • 10~*, we take B = 5, K = 15, 
a = 1, cf. Fig. 13.31 For n > 20 there is again no essential difference between the 
different Runge-Kutta methods. 

Fig. 13.41 shows the relative errors on the interval [10,20000] (similar for any in- 
terval [a, 2000a] with a > 10) for the RadauIIA(3) method with K = 10,20,40,80, 
160, 320, 640. For the imphcit Euler, the BDF(2) and the RadauIIA(5) method these 
error plots look similar. This behavior of the errors clearly demonstrates the ad- 
vantage of using local approximations. With B = 10, we need three approximation 
intervals to cover the interval [10, 20000], so that for a work oi 3 ■ K with if = 10 we 
obtain better accuracy than with K = 640 over the whole interval. 

In this example the maximum quadrature errors using the hyperbolas are smaller 
than those for the Talbot contours. Moreover, the hyperbolas allow to choose larger 
intervals. On the other hand, the Talbot contours turned out to be less sensitive to 




Fig. 3.4. Talbot (left) and hyperbola (right) quadrature error versus n for the RadauIIA(3) 
method with K = 10, 20, 40, 80, 160, 320, 640. The bold parts of the error curve correspond to the 
lower left parts of Figs. \3.'A and \3.S\ Note the different scaling of the y axis. 



the choice of parameters and the Laplace transform functions than the hyperbolas. 

3.3. Theoretical error bounds of the contour integral approximations. 

For the case of the hyperbola, we obtain in the same way as in Theorem 3 of the 
following error bound which shows exponential convergence. 

Theorem 3.1. There are positive constants C, d, cq, . . . ,C4, and c such that at 
t = nh < T the quadrature error in (15*.^ for a hyperbola ^3.51 is bounded by 

C COflt 



n/2 

if n > c^t and /it > 1 . Here v is the exponent of ^1.4[ ■ 

Given an error tolerance e, the first term in the error bound becomes 0{e ht'^~^) 
if T is chosen so small that co/it — iird/r < loge, which requires an asymptotic 
proportionality i ~ log - + fit. For /i chosen such that ^ log - < f^t < ai log - with 
an arbitrary positive constant oi and with _B > 1, we obtain that the second term is 
0{e ht'^~^) if ci — C2 cosh(i4'T) < —B/ai, i.e., with cosh(i^T) = 02 for a sufficiently 
large constant 02- With the above choice of r, this yields K ~ log K The third term 
then becomes smaller than e /it""^ for n> c log - with a sufhciently large constant c. 
In summary, this gives the following bound for the required number of quadrature 
points on the hyperbola. 

Theorem 3.2. In (j.'i'.ijjl . a quadrature error hounded in norm by eht'^~^ for 
nh e h is obtained with K — O(logi). This holds for n > c log i (with some 
constant c > 0), with K independent of £ and of n and h with nh < T . 

The approximation is, however, poor for the first few n, as we have seen in the 
numerical experiments. 

We refer to jS] for an optimized strategy to choose the parameters /i, r, K, which 
takes also perturbations in the evaluations of the Laplace transform into account. 

We expect that a similar result to Thcorcm l3 . 21 holds also for the Talbot contours, 
if the Laplace transform has an analytic continuation beyond the negative real axis 
from above and below, as is the case for the fractional powers considered above. 



4. The fast and oblivious algorithm. We now describe the convolution algo- 
rithm, concentrating on Runge-Kutta based convolution quadrature. The algorithm 
differs slightly depending on whether we want to compute a convolution or to solve 
an integral or integro-difFerential equation of convolution type. 

4.1. The algorithm for computing convolutions. The algorithm presented 
here uses the organisation scheme of the fast convolution algorithm described in a 
step by step manner in 03! ■ ^ pseudo-code for the algorithm developed in can 
be found in [S]. 

For fixed integer n < N and a given base B we split the discrete convolution H1.2(l 
or H2.7|l into L sums, where L is the smallest integer such that n < 2B^: 

n 

E(0) I , (L) 

^n-j gj = Kii + ■■■ + K+i 

bt-i-1 

with u^^l^ = ujo 9n and = ^ u/n-j 9j 

j=bt 

for suitable n = bo > bi > ■ ■ ■ > bh-i > 6l = 0. In view of the approximation 
intervals the splitting is done in such a way that for fixed I in each sum from 

bt to bi-i - 1, we have n - j e [S^^^, 2B^ - 2]. The bg = b(^), £ 1, . . . , L - 1 are 
determined recursively by the following pseudo-code. 

L = 1; q = 0; 
for n = 1 to N do 

if 2*B~L == n+1 then L = L+1; endif 

k = 1; 

while mod(n+l,B"k) == & k < L 

q(k) = q(k)+l; k = k+1; 
endwhile 

for k = 1 to L-1 do b(k) = q(k)*B~k; endfor 
endf or 

Note that for growing n, bg is augmented by every B^ steps. On inserting the 
integral representation (|2.9|l of the Runge-Kutta quadrature weights and the rela- 
tion (|2.1Q|I . i.e., e„_j(/iA) = r{hX)"^^^ q{hX), we obtain 

= J2 ^^-J9j= Y~- ^n-Ah\)®F{\)d\gj (4.1) 

= — f r(;iA)"-(''^-i-i'F(A)y(^H^A)dA 

with 

bf-l-l 

yW(/iA) = /i e(fc,_i-i)-j(/iA)<7j. 

j=bt 

Comparing this formula with H2.11|l . we see that y'^^\hX) is the Runge-Kutta approx- 
imation to the solution at i = bf^ih of the linear initial- value problem 
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(4.2) 



and hence y'^^\h\) is computed as such, by Runge-Kutta time-stepping. The integrals 
are discretized with the quadrature formula discussed in Section 13 

E w'i\{hXfr-''-^^^F{x'^^^)y^^\hX^^\ (4.3) 

k=-K 

In the nth time step, we thus compute and for subsequent use we update the 

Runge-Kutta solutions to the {2K+1)L initial value problems 14.2|l for the integration 
points Aj, on the contours for i = 1, . . . , L, doing one time step from t„ to tn+i 
in each of these differential equations. 

This algorithm does not keep the history gj {j = 0, . . . ,n) in memory. For each 
£ = 1, . . . ,L and k — —K, . . . ,K, it stores the Runge-Kutta approximation to H4.2|l 
at the current time step, the values w^i^\ x'f\ r{hxf^), F{X^^^), y^^'^ {hX^^^)^ and two 
auxiliary values of the dimension of y needed for book-keeping purposes (cf. |14[ I5]1. 
There are only {2K + 1)L evaluations of the Laplace transform F{s). In the case of 
real functions f{t) and g{t) only the real parts of the above sums are needed, and 
hence the factor 2K + 1 can be replaced by K + 1, since the quadrature points lie 
symmetric with respect to the real axis. We recall L < log^ N and K = 0(log -), 
where e is the accuracy requirement in the discretization of the contour integrals. 

In view of the poor approximation of the first convolution quadrature weights by 
the discretization of the contour integral, we evaluate u^^^i directly for a few of the 
first £, e.g., for £ ^ 0,1 with B = 10. For this we need to keep the n — 6i -|- 1 < 2B 
values gb^ , . . . , g„ in memory, but none of the earlier history gj for j < n — 2B. We also 
need the few convolution quadrature weights ujq, . . . , uj2B-i, which may be computed 
from (|3.3|l with 2B evaluations of the Laplace transform F{s). 

For the convolution quadrature based on the second-order BDF method a similar 
fast algorithm is obtained by inserting the formula (|2.4(l for e„-j(/iA) in H4.1|l . 

4.2. The algorithm for solving integral equations. The adaptation of the 
above algorithm to integral equations such as 

u{t)^ait)+ f fit-T)g{T,u{T))dT, t>0, (4.4) 
Jo 

is straightforward for the case of the convolution quadrature based on the implicit 
Euler method and the second-order BDF method, which use solution approximations 
only on the grid t = nh. The extension of the Runge-Kutta based algorithm is, 
however, less immediate, because the integral approximation uses the internal stages 
of the Runge-Kutta method. Consider a Runge-Kutta based convolution quadrature 
under the assumptions of Section 12.21 With the column vector of internal stages 
Vn — {vni)iLi, the discrctization of (|4.4|l reads 

n 

Vn = an + ^ Wn-j 5j , n > 0, (4.5) 

j=o 

with a„ = {a{tn + Ci/i))™ with weight matrices Wn defined by (|2.ti|) . and with 

9] ~ {di^j ~^ depending on the stages Vji. The scheme is implicit in Vn- 

The solution at i„+i is approximated by the last component of the stage vector Vn, 



^n+l Vnni 
10 



With the proof of 13, Theorem 4.1] we obtain that the error of this approximation 
over bounded time intervals is bounded by 0{h'^) with k — min(p, q+ 1), where p and 
q are the classical order and stage order, respectively, of the underlying Runge-Kutta 
method. This estimate holds under the assumption that the solution is sufficiently 
smooth. It gives orders 3 and 4 for the 2- and 3-stage Radau IIA methods, respectively. 
The precise approximation order for the 3-stage method (of classical order 5) may 
become larger under appropriate conditions on the nonlinearity and the convolution 
kernel, cf. [21 Theorem 4.2]. 

The weight matrix Wn has the integral representation, cf. (|2.9(l . 



Wn^— E^{h\) ® F 
2in Jp 



(A) dX , 

where the m x m matrix En{z) is defined by (|2.8|l . By Lemma 2.4 of ^Sli for n, > 1, 
En{z) is the rank-1 matrix given by 

En{z) = r(z)"-i(/ - za)-Hh^{I - zO)-^ 
= r(z)-i(/-za)-ille„(z). 

These relations permit us to proceed for the history term of l|4.5(l as we did for (|2.7|l . 
We split the stage vector w„ as 



with v\^^ = ^ Wn-j : 



-t;l"^H hwi^^ witn wr' = > w,,,-, gj 



and obtain, like in H4.1|l . 

= / rihX}''-'"-^ (I - h\a)-H(g) F(X) y^'^'HhX) dX , 

where y^^\hX) is again the Runge-Kutta approximation at t — bi^ih to the initial- 
value problem H4.2|l . now for the inhomogeneity values gj — {g{tj + Cih,Vji)^ in 

place of gj = {g{tj + Cih)Y^^. For ^ > 2 or 3, we thus approximate Vn^ as 

-i''- E 4^^(/.A[^^r-^-M/-^A[^^a)-^ii«^^(Ai^^)2/^^H^Ai^^). 

k=-K 

The algorithm stores the same values as before. The memory requirements for the 
algorithm are thus independent of the number of stages m and remain essentially the 
same as in the pure convolution case. 

5. Numerical experiments. We give two examples to illustrate the application 
and behavior of the fast convolution algorithm. 

5.1. A nonlinear Volterra equation. We consider a nonlinear Volterra inte- 
gral equation with weakly singular kernel from 

(z.(r)-sin(r))\ 

u{t)^- . , , — dr. (5.1) 

Jo y'n[t-T) 
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The convolution quadrature based on the backward Euler method gives the imphcit 
discretization 



Un = ^ujn-]{uj - sin(j/i))^, 
j=o 

where w„ is given by (|1.3|l with F{s) = s~^/^ and (5(C) — 1 — (. To solve the nonlinear 
equation in each time step we use Newton iterations. The history term is computed 
by the fast algorithm of the previous section. 

We consider also the discretizations based on the backward differentiation method 
of order 2, cf. Section [TTl and on the 2- and 3-stage RadauIIA implicit Rungc-Kutta 
methods of orders 3 and 5, respectively; see Sections 12.21 and |4 . 21 

In the numerical experiment we use the base B = b and the Talbot contours 
with K = 15 and K = 30 and the further parameters as in Section IX^ We choose 
a tolerance of 10^^^ in the Newton method. The error is calculated with respect to 
a reference solution, obtained with h — 0.001. Figure |0] shows the evolution of the 
absolute error and the oscillating solution u. 

Figure shows the errors u„ — u(t„) versus the step size h at time tn = 60, for 
K = 15 a,iidK ^ 30. 

Figure 1^31 plots the cpu time versus the number of integration steps, up to 10^ 
time steps. The near-linear growth of the computational work is clearly visible. 
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Fig. 5.1. Evolution of the solution over the interval [0,60] (left) and absolute error for different 
time integration methods (right) for h = 0.05 and K = 30. 




5.2. Fractional diffusion with transparent boundary conditions. Here 
we consider a fractional diffusion equation on the real line; see, e.g., jl5) for applica- 
tions of such equations in physics and for numerous references. The equation can be 
formulated as 

/■* (t - t)""1 

u(x,t) - uq(x) = / dxxu(x,T)dT + g(x,t) {otxgR, t>0 (5.2) 

Jo r(a) 

with the asymptotic condition u{x, i) — > for x — > ±cx), for an inhomogeneity g with 
gi^XjO) = 0. To reduce the computation to a finite domain x € [—a, a] for initial data 
Wo and inhomogeneity g with support in [—a, a], we impose transparent boundary 
conditions at x = ±a, which read 

u{±a,t) = ~ L_^^_a^^,(±a,^)d^, (5.3) 
Ja r(a/2) 
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Fig. 5.2. Absolute error vs. step size h, for different integration methods, with K = 15 (left) 
and K = 30 (right). 
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Fig. 5.3. Cpu time in seconds versus the number of integration steps. 



with the outward derivative di, — ±dx at x = ±a. These boundary conditions are 
derived with Laplace transform techniques in the same way as for the wave or the 
Schrodinger equation; see, e.g., |[2j. Space discretisation of H5.2(l is done using sec- 
ond order finite differences and a central finite difference to approximate the normal 
derivative. With the notation 

^xxu" = Aa!^^"""^ ^ ^ """+1) ' '^i^"±(M-i) — QAx^^^^^^ ^ ^±{M-2)) 
for a = MAx, the discrete equation approximating H5.2|) is 

n 

= 4-1 a" for( = -(A/-l),...,Af-l ; n > 0, 

3=0 
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where the weights ojfi are the convohition quadrature weights for the kernel f{t) = 
t^~^/r{(3) with Laplace transform F{s) = s^^. 

In the numerical example we set a ~ 5 and A/ = 450. We consider the problem 
with a — 2/3 and no inhomogeneity, i.e., g = 0. The initial value is u{x,0) = 
exp(— x^). Figure 1^31 shows the errors at t = 2 in dependence on the step size for 
the Radau IIA methods of orders 1, 3, 5, obtained with B = 5 and K = 15 in the 
fast convolution algorithm. The reference solution is obtained with the Radau IIA 
method of order 5, with h = 0.0002 and K = 40. We observe an order reduction 
for the higher-order methods, which is due to the temporal non-smoothness of the 
solution at i = 0; cf. jlj Sect. 8]. Nevertheless, the higher-order methods give much 
better accuracy. 

The work diagram looks almost identical to Figure [^31 showing practically linear 
dependence of the computational work on the number of time steps. The required 
memory is less than 200 entries per spatial grid point for up to A'^ < 10** steps, and 
less than 300 entries per grid point for A^ < 10^ steps. These numbers are halved if 
we run the algorithm with B = 10, K = 10 instead oi B = 5, K = 15, as is sufficient 
for less stringent accuracy requirements (~ 10~^). 




step size 

Fig. 5.4. Absolute error vs. time step, for different integration methods, with K = 15. 
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